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Generation of random thermal particle momenta is a basic task in many problems, such as microscopic studies of 
equilibrium and transport properties of systems, or the conversion of a fluid to particles. In heavy-ion physics, the 
(in)cfficicncy of the algorithm matters particularly in hybrid hydrodynamics + hadronic transport calculations. With 
popular software packages, such as UrQMD 3.3pl [1] or THERMINATOR @], it can still take ten hours to generate 
particles for a single Pb+Pb "event" at the LHC from fluid dynamics output. Below I describe reasonably efficient 
basic algorithms using the MPC package (Molnar's Program Collection) Q, which should help speed momentum 
^—^ ' generation up by at least one order of magnitude. 

. It is likely that this wheel has been reinvented many times instead of reuse, so there may very well exist older 
and/or better algorithms that I am not aware of (MPC has been around only since 2000). The main goal here is to 

\ encourage practitioners to use available efficient routines, and offer a few practical solutions. 
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ON ■ A. Static thermal distributions 

^ c"| First consider a static thermal distribution 
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where /x is the chemical potential, m is the particle mass, g is the degeneracy factor (internal degrees of freedom), and 
I . C +1 for fermions 

■ a = < — 1 for bosons (2) 
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and the only challenge is to generate the scaled magnitude x — \p\/T 

. . dN x 2 m /T 

A straightforward technique to employ is an automated version of the rejection method (see Fig. Q}. The rejection 
method 4] is based on a suitable comparison function F that bounds g from above for all x. One then generates 
uniformly (x, y) pairs in the two-dimensional area under the graph of F, and keeps the x coordinate of those points 
that are below the graph of g. It is easy to construct a staircase comparison function, automatically, via tabulating 
g(x) at several points, and multiplying by a modest "safety" factor to leave room for any local maxima within the 
intervals Uniform sampling under the graph of a staircase function is straightforward Q 

Alternatively, we can tabulate g(x), construct an interpolation of it, and generate the interpolated distribution [§[. 
Computationally this is very similar to staircase function sampling in automated rejection. To precisely represent 
g{x), in general many points are needed but we then have 100% efficiency (there is no rejection check). 



B. Boosted thermal distributions 



For thermal distribution with collective flow 



/ - (p) '^ = (2^ e(^-^ + a ' P° = ^^ • «"=7(l,v,) (5) 

where v> is the three-velocity of the fluid. This is a special case (fi/ 1 = (1,0)) of Cooper-Frye freezeout discussed in 
the next Section. 
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FIG. 1: Illustration of rejection method with staircase comparison function. 



C. Cooper- Prye freezeout 



In Cooper-Frye freezeout [6[ a fluid is converted to particles across a three-dimensional spacetime hypersurface. A 
surface element do~^(x) at spacetime point x contributes to particle number 



dN =p a da a (x)^2f(x,p) 
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where / is the phase space density of the particles in the fluid, 
distribution of particles emitted from the surface element is 

dN p a n c 



For an ideal fluid / is thermal, so the momentum 



fv(p) 



with /„ from Because dN is a Lorentz scalar, in the local rest frame of the fluid element (v> 
static thermal distribution with a momentum-dependent prefactor 
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After generating a thermal momentum p, the prefactor can be accounted for in an additional rejection step. If p fails 
the rejection test, we generate a new one, otherwise we boost it back to the original frame to obtain p. 

For timelike normal with n° > 0, the prefactor is always nonnegative, and the additional rejection step is at least 
50% efficient because we can use as comparison function h° + |n| < 2h°, while on average (n° — nv) = n°. Negative 



< is unphysical (dN is never positive then). So we need no more than two tries, on average. 



For spacelike normal (n 2 = — 1), there are always momenta for which (n° — nv) is positive, and also momenta for 
which it is negative, so one needs to follow some prescription to handle unphysical contributions. A customary choice 
is to ignore negative values and implement a cut prefactor (n° — nv)0(n° — nv) instead, even if with thermal fo 
conservation laws are violated^. For h° > 0, a comparison function h° + |n| then provides II) an efficiency of at 
least (|v|)/4, which for ra/T ^ 1 is poor but still workable in heavy-ion physics applications because spacelike surface 
elements are rare, and heavy particles have low abundances. For n° < the total particle yield from the hypersurface 
element (without the cut) is negative, and for that reason such surface elements are most often omitted. Momentum 
generation for n° < 0, if desired, requires additional care fill]. 

For better efficiency, especially for spacelike normals and heavy particles, we can generate directly the isotropic 
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and realize the correct prefactor (n° — nv)0(n° — nv) via rejection using the more efficient |12| comparison function 
n° + |n||v|. This way, for timelike normals we need no more than two tries, on average, while for spacelike normals 
with n° > 0, no more than four tries. The price is that a new generator has to be constructed for each hypersurface 
element, whereas the static thermal generator (/o) can be reused when T = const across the entire hypersurface (often 
the case in practice). 
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D. Massless Boltzmann distribution 



Finally, I show that for massless Boltzmann particles static thermal momenta 

f **e- , (10) 
can be generated through the one-liner [l3j 

^ = -ln[(l-ei)(l-6)(l-6)] (11) 

where £1,2,3 are three independent uniform deviates on [0,1). This follows because with x — r 2 , the distribution 
corresponds to a Gaussian in six-dimensional space {y\, yz, . . . , ye} 

dN 



r 5 dr 



(xe- r , r = Jy( + yi + --- + yi (12) 



which can be generated as three independent pairs of Gaussian deviates, each pair through the standard 2D polar 
coordinate technique dx dy e~ x e~ y oc dip d(p 2 ) e~ p . Here ip is uniform on [0, 2ir) and p 2 is exponential on [0, 00). We 
only need p 2 for each pair because r 2 = p\ 2 + p| 4 + p 2 6 , and exponential deviates can be generated through inverting 



I{p ) = j dwe' w = l-e- p =S> p z = -\n(l-I) (13) 



where / is a uniform deviate on [0, 1). 



E. Conclusion 

In a paradigm where random numbers are computationally very expensive, a useful efficiency benchmark for three- 
momentum generators is 

(14) 



(uniform deviates per momentum) 

(i.e., with one random number generator call per momentum component, e = 1). For the algorithms above - static 
thermal (Sec. A), boosted thermal (Sec. B), and Cooper-Frye (Sec. C) - e = 1, e > 3/8, and e > 3/16, respectively, 



if we use the interpolation method [14 1 



In reality there is of course overhead due to function tabulation, memory lookups, and binary search, which limit 
the number of intervals one can utilize. But in practice, overall efficiency is still a reasonable ~ 0.1 (even with the 
automated rejection method), and often much better [l5|. 
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